% code for robust averaging of SPMs
% Ref:1) Main robust- Ritter, O., A. and Dawes, G.J.K., 1998. New equipment and processing 
%        for magnetotelluric remote refernce observations, GJI. 132, 535-548.   
% 		2) Jones et al 1989, A comparison of techniques for magnetotelluric 
%        response function estimation, JGR, 94, 14201-14213.
%		3) Bivariate Error- Muller, A. 2000. A new method o compensate for bias
% 			in magnetotellurics Geophy. J. Int 142, 257-269
% Tested OK 12.4.3
%
% development history
% 20.3.3 coding according to Oliver Ritter, 1998.
% 5.4.3 initial tf is now computed from averaged SPM
% 22.4.3 initial tf is computed from medain SPM
% 23.4.5 Tukey weights added to reduce extreme outliers (Beaton & Tukey, 1974)
% 24/4  deleted the second weighting step (Lc etc)% reinstated 1.5.2003
% 25.4.3 variance estimations acc to Muller, 2000
% 10.6.3 Jack Kniving as inital guess
% Latest date 30.06.2003

function[tf_t,ProcDef,SP] = robspm(SPMall,ProcDef),

fprintf('Robspm-> Robust estimation of transfer functions\n');
A = size(SPMall);
ProcDef.Weights=ones([A(1),A(2)]);

if ProcDef.overlap ==1,  % Overlap of windows cause reduction in dof for FC's
	v=0.818;					 % A correction to this effect is done here
else,
 	v=1;
end;

%SPMall = NormSPM(SPMall); 		% Normalizing SPM wrt power in horizontal magnetic fileds
fprintf('Jack kniving for initial guess...\n');
[z,dof]=jknf(SPMall,ProcDef);

fff = A(1)/12;
aaa = 1;
fprintf('0         100percent\n');

for i = 1:A(1),
      
   if i/fff > aaa,
     		fprintf('#');
      	aaa = aaa+1;
   end;

   
  SPMat(:,:,:) = SPMall(i,:,:,:); % all the events for one frequency
   
% Inserted on 13.06.2003 Here a number of algorithms for initial guess can be thought of
% LS, median, coherence & Jack knife are currently used ( see Egbert & Livelybrooks Geophysics ,Chave et al 1989)
%tfmed = tf1(SPMat,ProcDef);			% All event tfs, single frequency 
%ti = median(tfmed); 					% single tf matrix (median) (initial guess)
%SPMi=cohguess(SPMat,ProcDef,1,i);
%ti=tf(SPMi,ProcDef);
ti= z(i,1:4);  
for ii = 1:A(2),
	data(:,:) = SPMat(ii,:,:);  		% s
   a = abs(r_resid(ti,data,'Ex'));
   b = abs(r_resid(ti,data,'Ey'));
   S(ii) = sqrt(a^2+b^2);            % residual to median tf
   S(ii)=a;
end;											% the residual returned is |EobsEobs* - EpreEpre*|
	Temp=S;
   sig_m = 1.483*median(abs(S-median(S))); % median absolute deviation 
   c_m = 1.5*sig_m;
  
   k = find(S <= c_m);					% assigning Huber weights 
   w_m(k) = 1;
   q = find(S > c_m);
  	w_m(q) = c_m./S(q);

   ProcDef.Weights(i,:) = w_m;
   
   dof_factor = sum(w_m)/A(2); % effect of weighting on degrees of freedom
   
	[SPMi,SPMat] = WeightSPM(SPMat,ProcDef,i); % weigthing matrices

   tf_m  = tf1(SPMi,ProcDef);
   
   for ii = 1:A(2),
   	data(:,:) = SPMat(ii,:,:);
      a = abs(r_resid(tf_m,data,'Ex'));
      Temp1(ii)=r_resid(tf_m,data,'Ex');
	   b = abs(r_resid(tf_m,data,'Ey'));
      S(ii) = sqrt(a^2+b^2);
      S(ii)=a;
   end;
   
   %-----------------MORE SEVERE WEIGHTS---------------------------
   Lc = length(k);
   
   if Lc < 1,
      fprintf('robspm-> Warning : Adjust c_m value\n');
      Lc = 1;
   end;
   
  	sig_h = (A(2)/Lc^2)*sum((w_m.*S).^2); % new variance estimate
   c_h = 0.5*sqrt(sig_h);% Sqrt added on 7.4.3 as in A1.3 Ritter et al

   
   k = find(S <= c_h);		% new weights
   w_h(k) = 1;
   q = find(S > c_h);
   w_h(q) = c_h./S(q);
 
   ProcDef.Weights(i,:) = w_h;
   dof_factor = sum(w_h)/dof_factor; % re weighting & degrees of freedom
   [SPMi,SPMat] = WeightSPM(SPMat,ProcDef,i); % weigthing matrices
   tf_h(i,:)  = tf1(SPMi,ProcDef);
   
   for ii = 1:A(2),
	   data(:,:) = SPMat(ii,:,:);
   	a = abs(r_resid(ti,data,'Ex'));
   	b = abs(r_resid(ti,data,'Ey'));
      S(ii) = sqrt(a^2+b^2);            % Weighted SPM residual to new tf
      S(ii)=a;
      
   end;

   %----------------NOW GET RID OF EXTREME OUTLIERS BY TUKEY WIEGHTS---------------
   
   sig_t = abs(sqrt((1/A(2))*sum((w_h.*S).^2)/...
   ((1/Lc)*sum((1-(S./c_h).^2).*(1-5*(S./c_h).^2))))); % Junge
 	c_t = 6*sig_t;
   
   S=Temp; % inserted 26.4.3 by passes the stage 2
   k = find(S <= c_t);		% Tukey's biweights 
   w_t(k) = (1-(S(k)./c_t).^2).^2;
   q = find(S > c_t);
   w_t(q) = 0;
   ProcDef.Weights(i,:) = w_t;
   
   dof_factor = sum(w_t)/dof_factor;
   
   [SPMi,SPMat] = WeightSPM(SPMat,ProcDef,i); % weigthing matrices
	tf_t(i,1:4)  = tf1(SPMi,ProcDef);
   SP.spectra(i).data=squeeze(SPMi); %squeeze added on 18.8.04
   SP.spectra(i).freq = ProcDef.TLFreq1(i);
%------------------Estimate the BIVARIATE ERROR------------------------------
	

	data(:,:)=SPMi(1,:,:);
	c1 = cohe(tf_t(i,3),tf_t(i,1),data,'Ex'); % predictive coherency ex-hxhy
	c2 = cohe(tf_t(i,2),tf_t(i,4),data,'Ey');% predictive coherency ey-hxhy
   dof = ProcDef.TLRad1(i)*v*2*dof_factor*A(2); % degrees of freedom = frequency band width*nsegments*weight factor
   SP.spectra(i).dof.avgt=dof;
   SP.spectra(i).dof.avgf=1;
	[tf_t(i,7),tf_t(i,5)] = BivError(c1,dof,data,'Ex'); % bivaraite error
	[tf_t(i,6),tf_t(i,8)] = BivError(c2,dof,data,'Ey');

end;
SP.head.site = ProcDef.file1;
fprintf('\n Over..Switching to main\n');

